{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 习题 9.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如例9.1的三硬币模型，假设观测数据不变，试选择不同的处置，例如，$\\pi^{(0)}=0.46,p^{(0)}=0.55,q^{(0)}=0.67$，求模型参数为$\\theta=(\\pi,p,q)$的极大似然估计。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class EM:\n",
    "    def __init__(self, prob):\n",
    "        self.pro_A, self.pro_B, self.pro_C = prob\n",
    "        \n",
    "    def pmf(self, i):\n",
    "        pro_1 = self.pro_A * math.pow(self.pro_B, data[i]) * math.pow((1-self.pro_B), 1-data[i])\n",
    "        pro_2 = (1 - self.pro_A) * math.pow(self.pro_C, data[i]) * math.pow((1-self.pro_C), 1-data[i])\n",
    "        return pro_1 / (pro_1 + pro_2)\n",
    "    \n",
    "    def fit(self, data):\n",
    "        print('init prob:{}, {}, {}'.format(self.pro_A, self.pro_B, self.pro_C))\n",
    "        count = len(data)\n",
    "        theta = 1\n",
    "        d = 0\n",
    "        while(theta > 0.00001):\n",
    "            # 迭代阻塞\n",
    "            _pmf = [self.pmf(k) for k in range(count)]\n",
    "            pro_A = 1/ count * sum(_pmf)\n",
    "            pro_B = sum([_pmf[k]*data[k] for k in range(count)]) / sum([_pmf[k] for k in range(count)])\n",
    "            pro_C = sum([(1-_pmf[k])*data[k] for k in range(count)]) / sum([(1-_pmf[k]) for k in range(count)])\n",
    "            d += 1\n",
    "            print('{}  pro_a:{:.4f}, pro_b:{:.4f}, pro_c:{:.4f}'.format(d, pro_A, pro_B, pro_C))\n",
    "            theta = abs(self.pro_A-pro_A) + abs(self.pro_B-pro_B) + abs(self.pro_C-pro_C)\n",
    "            self.pro_A = pro_A\n",
    "            self.pro_B = pro_B\n",
    "            self.pro_C = pro_C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载数据\n",
    "data=[1,1,0,1,0,0,1,0,1,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "init prob:0.46, 0.55, 0.67\n",
      "1  pro_a:0.4619, pro_b:0.5346, pro_c:0.6561\n",
      "2  pro_a:0.4619, pro_b:0.5346, pro_c:0.6561\n"
     ]
    }
   ],
   "source": [
    "em = EM(prob=[0.46, 0.55, 0.67])\n",
    "f = em.fit(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可见通过两次迭代，参数已经收敛，三个硬币的概率分别为0.4619，0.5346，0.6561"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 习题9.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "已知观测数据  \n",
    "-67，-48，6，8，14，16，23，24，28，29，41，49，56，60，75  \n",
    "试估计两个分量的高斯混合模型的5个参数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## sklearn.mixture.GaussianMixture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__参数：__ \n",
    "- __n_components：__ 混合高斯模型个数，默认为1 \n",
    "- __covariance_type：__ 协方差类型，包括{‘full’,‘tied’, ‘diag’, ‘spherical’}四种，分别对应完全协方差矩阵（元素都不为零），相同的完全协方差矩阵（HMM会用到），对角协方差矩阵（非对角为零，对角不为零），球面协方差矩阵（非对角为零，对角完全相同，球面特性），默认‘full’ 完全协方差矩阵 \n",
    "- __tol：__ EM迭代停止阈值，默认为1e-3. \n",
    "- __reg_covar：__ 协方差对角非负正则化，保证协方差矩阵均为正，默认为0 \n",
    "- __max_iter：__ 最大迭代次数，默认100 \n",
    "- __n_init：__ 初始化次数，用于产生最佳初始参数，默认为1 \n",
    "- __init_params：__ {‘kmeans’, ‘random’}, defaults to ‘kmeans’.初始化参数实现方式，默认用kmeans实现，也可以选择随机产生 \n",
    "- __weights_init：__ 各组成模型的先验权重，可以自己设，默认按照init_params产生 \n",
    "- __means_init：__ 初始化均值，默认按照init_params产生\n",
    "- __precisions_init：__ 初始化精确度（模型个数，特征个数），默认按照init_params实现 \n",
    "- __random_state：__ 随机数发生器 \n",
    "- __warm_start：__ 若为True，则fit（）调用会以上一次fit（）的结果作为初始化参数，适合相同问题多次fit的情况，能加速收敛，默认为False。 \n",
    "- __verbose：__ 展示迭代信息，默认为0，可以为1或者大于1（显示的信息不同） \n",
    "- __verbose_interval：__ 若展示迭代信息，设置多少次迭代后显示信息，默认10次。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.mixture import GaussianMixture\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 初始化观测数据\n",
    "data=np.array([-67, -48, 6, 8, 14, 16, 23, 24, 28, 29, 41, 49, 56, 60, 75]).reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "labels = [0 0 1 1 1 1 1 1 1 1 1 1 1 1 1]\n"
     ]
    }
   ],
   "source": [
    "# 聚类\n",
    "gmmModel = GaussianMixture(n_components=2)\n",
    "gmmModel.fit(data)\n",
    "labels = gmmModel.predict(data)\n",
    "print(\"labels =\", labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGeNJREFUeJzt3Xm4XXV97/H3R46KQDBgQhiCBhUHijTSIxf0lmsFKiKCz72mxapQ9ZY6a9UroPdxemqd+mj12lojqFCplOBA7NUqopbbimjQiAIqKVQIQzgOQJzAA9/7x1pHNuFknUOSnXWG9+t58uy9hr3Xd59zsj/r9/utIVWFJEmbc7++C5AkzWwGhSSpk0EhSepkUEiSOhkUkqROBoUkqZNBoTkpyeeTnNTj9h+a5OdJduirhu0tyZOTrJ/mum9O8vFh16Rtw6DQNpHkhCSXJPlFkpvb5y9Jkj7qqaqnVdWZ2/p9k/xpkkrynk3mP7Od/7F2+9dW1S5Vdec03vNjSf5yW9c6xTYryYYkIwPzRtrfnSdX6R4MCm21JK8B3ge8G9gTWAK8CHgS8IAeSxuW/wD+ePBLFjgR+GEfxWxSx31xC/C0geljgJ9tfUWaawwKbZUkDwbeCrykqs6rqo3V+HZVPaeqbm/Xe3qSbye5Lcl1Sd488B736rJI8p9JjmyfH5JkTfvaDRN780l2TPLxJD9JckuSbyZZ0i77apL/2T5/RJIvt+v9OMnZSRZusq3XJrksya1J/inJjh0f+ybgu8BT29fvDjwRWD3wnsvavfaRJLsnWZ/kGe2yXZKsS3JikpOB5wCva7uqPtuuU0keOfB+v211TPy8kpyS5Cbgo+38Y5OsbX8WX0ty0BS/vn+gCbgJJwJnbfJ72DvJ6iQ/bWv+s4FlD2rr+lmSK4AnTPLaTyYZS3JNkldMUY9mKINCW+sw4IHA+VOs9wuaL6KFwNOBFyd55jS38T7gfVW1K/AI4Nx2/knAg4F9gYfQtGJ+NcnrA7wd2Bt4bLv+mzdZ54+Ao4H9gIOAP52iprO4+0v2BJrPf/tkK1bVT4EXAB9OsgfwXmBtVZ1VVSuBs4F3tV1Vz5hiuxP2BHYHHgacnORg4CPAn9P8LD4ErE7ywI73+AxweJKFbXD+Pvf+PX4CWE/zs3sW8FdJjmiXvYnm9/EImtD87ZhQkvsBnwW+A+wDHAG8KslTp/n5NIMYFNpai4AfV9X4xIx2b/aWJL9KcjhAVX21qr5bVXdV1WU0X0D/bZrb+A3wyCSLqurnVfX1gfkPAR5ZVXdW1aVVddumL66qdVV1QVXdXlVjwHsm2fb7q+qG9kv9s8DyKWr6NPDktkV1rz3xSWr4IrAKuJAmKP98ivefyl3Am9rP9Cvgz4APVdUl7c/iTJrgOrTjPX5N81n/mCbsVrfzAEiyL/BfgVOq6tdVtRY4HXheu8ofAW+rqp9W1XXA+wfe+wnA4qp6a1XdUVVXAx9ut6NZxqDQ1voJsGiwn7yqnlhVC9tl9wNI8l+SfKXthriVZu9/0TS38ULgUcD32+6lY9v5/wB8ATgnyQ1J3pXk/pu+OMkeSc5Jcn2S24CPT7Ltmwae/xLYpaug9sv5/wL/G1hUVf8+jc+xEjgQ+GhV/WQa63cZq6pfD0w/DHhNG9C3JLmFpuW09xTvM9Eymizs9gZ+WlUbB+b9iKaFMLH8uk2WDdaz9yb1vJ5m/EqzjEGhrXUxzZ7r8VOs9480e6z7VtWDgb+n6RKCpltqp4kV0xxSunhiuqquqqpnA3sA7wTOS7JzVf2mqt5SVQfQjBEcyz373Ce8HSjgoLb76rkD294aZwGvoQmsTu1n+lD7mhcPjj+0tW3qlwz8TGi6mgZt+prraPbuFw7826mqPjFFaf8P2IvmC/zfNll2A7B7kgUD8x4KXN8+v5EmjAaXDdZzzSb1LKiqY6aoRzOQQaGtUlW3AG8B/i7Js9qB2vslWQ7sPLDqApq9018nOQT4k4FlPwR2bAe870+zl/7bvvUkz02yuKruojlSB+DOJH+Q5HHtl/BtNF1Rkx2OugD4OXBLkn2A/7VNPjz8K3AU8H+mse7r28cXAH8NnJW7z7HYADx8k/XXAn+SZIckRzN1N92HgRe1Lbck2bn9eS7oelE19xl4BnBcbXLPgbY76WvA29sDBw6iad2d3a5yLnBakt2SLAVePvDybwC3tQPuD2o/x4FJ7jHgrdnBoNBWq6p3Aa8GXgfcTPPF9yHgFJovGoCXAG9NshF4I3cPSFNVt7bLT6fZW/0FzQDqhKOBy5P8nGZg+4S222VP4DyakLiS5ot7spO43gIcDNxK0130qa3+0E3dVVUXtuMam5Xk92h+Pie251W8k6ZFcGq7yhnAAW0XzWfaea+k+QK/heaoqM/QoarW0IxTfIDmENd1TD0gP/Hay6vq8s0sfjawjKZ18WmacZEL2mVvoeluugb4IgMtq/ZzPoNmrOca4Mc0v98HT6cmzSzxxkWSpC62KCRJnQwKSVIng0KS1MmgkCR12tKLic0oixYtqmXLlvVdhiTNKpdeeumPq2rxVOvNiaBYtmwZa9as6bsMSZpVkvxo6rV67npK8hdJLk/yvSSfaE/q2S/NvQyuSnMVz7l4mWpJmjV6C4r2DNlXAKNVdSCwA80Fw94JvLeq9qc5ceiFfdUoSep/MHsEeFB7QbmdaK4d8xSas20BzgSmeylqSdIQ9BYUVXU9zTVvrqUJiFuBS4FbBi5ZvZ67r1R5D0lOTnMzmzVjY2Pbo2RJmpf67HrajeaKo/vRXK54Z+55W8YJk15jpKpWVtVoVY0uXjzloL0kaQv12fV0JM1liMeq6jc0F2p7IrBw4N4GS2kuRiZJ6kmfQXEtcGiSnZKE5laJVwBfobnlIjS3VpzqFpuSpCHqc4ziEppB62/R3Kj+fjR3ADsFeHWSdTS3uTyjrxolST0f9VRVb6qqx1TVgVX1vPb+v1dX1SFV9ciqWlFVk96wXpLmow0b4PDDYdddm8cNG4a/zb4Pj5Uk3QcrVsDFF8PGjc3jihXD36ZBIUmzyNq1MN6eQDA+3kwPm0EhSbPI8uUw0h4XOjLSTA+bQSFJs8iqVXDYYbBgQfO4atXwtzknrh4rSfPFkiVw0UXbd5u2KCRJnQwKSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQpLUyaCQJHUyKCRJnQwKSVKnXoMiycIk5yX5fpIrkxyWZPckFyS5qn3crc8aJWm+67tF8T7gX6rqMcDvAlcCpwIXVtX+wIXttCTNKhs2wOGHw667No8bNvRd0ZbrLSiS7AocDpwBUFV3VNUtwPHAme1qZwLP7KdCSdpyK1bAxRfDxo3N44oVfVe05fpsUTwcGAM+muTbSU5PsjOwpKpuBGgf95jsxUlOTrImyZqxsbHtV7UkTcPatTA+3jwfH2+mZ6s+g2IEOBj4YFU9HvgF96GbqapWVtVoVY0uXrx4WDVK0hZZvhxG2ptNj4w007NVn0GxHlhfVZe00+fRBMeGJHsBtI8391SfJG2xVavgsMNgwYLmcdWqvivacr0FRVXdBFyX5NHtrCOAK4DVwEntvJOA83soT9I8MoyB5yVL4KKL4LbbmsclS7b+Pfsy0vP2Xw6cneQBwNXA82nC69wkLwSuBWbxEJCk2WBi4Hl8/O6B54su6ruqmaPXoKiqtcDoJIuO2N61SJq/5tLA8zD0fR6FJPVuLg08D4NBIWnem0sDz8PQ9xiFJPVuYuBZk7NFIUnqZFBIkjoZFJKkTgaFJKmTQSFp1phLl+6eTQwKSbPGXLp092xiUEiaNTyDuh8GhaRZwzOo+2FQSJo1PIO6H56ZLWnW8AzqftiikCR1MigkSZ0MCklSJ4NCktTJoJAkdeo9KJLskOTbSf65nd4vySVJrkryT+39tCXNMl5uY+7oPSiAVwJXDky/E3hvVe0P/Ax4YS9VSdoqXm5j7ug1KJIsBZ4OnN5OB3gKcF67ypnAM/upTtLW8HIbc0ffLYq/AV4H3NVOPwS4paraPy/WA/tM9sIkJydZk2TN2NjY8CuVdJ94uY25o7egSHIscHNVXTo4e5JVa7LXV9XKqhqtqtHFixcPpUZJW87LbcwdfV7C40nAcUmOAXYEdqVpYSxMMtK2KpYCN/RYo6Qt5OU25o7eWhRVdVpVLa2qZcAJwJer6jnAV4BntaudBJzfU4mSJPofo5jMKcCrk6yjGbM4o+d6JGlemxFXj62qrwJfbZ9fDRzSZz2SpLvNxBaFJGkGMSgkSZ0MCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQprnvG+EpmJQSPOc943QVAwKaZYY1p6/943QVAwKaZYY1p6/943QVAwKaZYY1p6/943QVGbERQElTW358qYlMT6+bff8vW+EpmKLQpol3PNXXwwKaQiGMfA8sed/223N45IlW/+e0nQYFNIQeMip5hKDQvOah5xKUzMoNK95yKk0NYNC85qHnEpT6y0okuyb5CtJrkxyeZJXtvN3T3JBkqvax936qlFz37D2/B141lzSZ4tiHHhNVT0WOBR4aZIDgFOBC6tqf+DCdloaCvf8pan1dsJdVd0I3Ng+35jkSmAf4Hjgye1qZwJfBU7poUTNA55sJk1tRoxRJFkGPB64BFjShshEmOyxmdecnGRNkjVjY2Pbq1RJmnd6D4okuwCfBF5VVbdN93VVtbKqRqtqdPHixcMrUJLmuV6DIsn9aULi7Kr6VDt7Q5K92uV7ATf3VZ8kqd+jngKcAVxZVe8ZWLQaOKl9fhJw/vauTTOTd2KT+tFni+JJwPOApyRZ2/47BngHcFSSq4Cj2mnJy2JIPenzqKd/A7KZxUdsz1o0O3hZDKkfvQ9ma+4ZVheRl8WQ+mFQaJsbVheRJ8dJ/fAOd9rmhtVF5MlxUj9sUWibs4tImlsMCm1zdhFJc4tdT9rm7CKS5hZbFJKkTgaFJKmTQTGPeUkMSdNhUMxjXhJD0nQYFPOYl8SQNB0GxSwxjG4iz3eQNB0GxSwxjG4iz3eQNB2eRzFLDKObyPMdJE3HlC2KJC9Lstv2KEabZzeRpL5Mp+tpT+CbSc5NcnR7ZzptZ3YTSepLqmrqlZpw+EPg+cAocC5wRlX9x3DLm57R0dFas2ZN32VI0qyS5NKqGp1qvWkNZleTJje1/8aB3YDzkrxrq6qUJM140xmjeEWSS4F3Af8OPK6qXgz8HvA/hlVY2831gyTrkpw6rO1IkrpN56inRcB/r6ofDc6sqruSHDuMopLsAPwtcBSwnmaMZHVVXTGM7UmSNm/KFkVVvXHTkBhYduW2LwmAQ4B1VXV1Vd0BnAMcP6RtSZI6zNQT7vYBrhuYXt/O+60kJydZk2TN2NjYdi1OkuaTmRoUkx2Ce4/Ds6pqZVWNVtXo4sWLt1NZkjT/zNSgWA/sOzC9FLihp1okaV6bqUHxTWD/JPsleQBwArC655okaV6akdd6qqrxJC8DvgDsAHykqi7vuSxJmpdmZFAAVNXngM/1XYckzXcztetJkjRDGBSSpE4GhSSpk0EhSepkUEiSOhkUkqROBoUkqZNBIUnqZFBIkjoZFJKkTgaFJKmTQSFJ6mRQSJI6GRSSpE4GhSSpk0EhSepkUEiSOhkUkqROBoUkqVMvQZHk3Um+n+SyJJ9OsnBg2WlJ1iX5QZKn9lGfJOlufbUoLgAOrKqDgB8CpwEkOQA4Afgd4Gjg75Ls0FONkiR6Coqq+mJVjbeTXweWts+PB86pqtur6hpgHXBIHzVKkhozYYziBcDn2+f7ANcNLFvfzruXJCcnWZNkzdjY2JBLlKT5a2RYb5zkS8Cekyx6Q1Wd367zBmAcOHviZZOsX5O9f1WtBFYCjI6OTrqOJGnrDS0oqurIruVJTgKOBY6oqokv+vXAvgOrLQVuGE6FkqTp6Ouop6OBU4DjquqXA4tWAyckeWCS/YD9gW/0UaMkqTG0FsUUPgA8ELggCcDXq+pFVXV5knOBK2i6pF5aVXf2VKMkiZ6Coqoe2bHsbcDbtmM5kqQOM+GoJ0nSDGZQSJI6GRSSpE4GhSSpk0EhSepkUEiSOhkUkqROBoUkqZNBIUnqZFBIkjoZFJKkTgaFJKmTQSFJ6mRQSJI6GRSSpE4GhSSpk0EhSepkUEiSOvUaFElem6SSLGqnk+T9SdYluSzJwX3WJ0nqMSiS7AscBVw7MPtpwP7tv5OBD/ZQmiRpQJ8tivcCrwNqYN7xwFnV+DqwMMlevVQnSQJ6CookxwHXV9V3Nlm0D3DdwPT6dp4kqScjw3rjJF8C9pxk0RuA1wN/ONnLJplXk8wjyck03VM89KEP3cIqJUlTGVpQVNWRk81P8jhgP+A7SQCWAt9KcghNC2LfgdWXAjds5v1XAisBRkdHJw0TSdLW2+5dT1X13arao6qWVdUymnA4uKpuAlYDJ7ZHPx0K3FpVN27vGiVJdxtai2ILfQ44BlgH/BJ4fr/lSJJ6D4q2VTHxvICX9leNJGlTnpktSepkUEiSOhkUkqROBoUkqZNBsa1t2ACHHw677to8btjQd0WStFUMim1txQq4+GLYuLF5XLGi74okaasYFNva2rUwPt48Hx9vpiVpFjMotrXly2GkPT1lZKSZlqRZzKDY1latgsMOgwULmsdVq/quSJK2Su9nZs85S5bARRf1XYUkbTO2KCRJnQwKSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktTJoJAkdTIoJEmdeguKJC9P8oMklyd518D805Ksa5c9ta/6JEmNXi7hkeQPgOOBg6rq9iR7tPMPAE4AfgfYG/hSkkdV1Z191ClJ6q9F8WLgHVV1O0BV3dzOPx44p6pur6prgHXAIT3VKEmiv6B4FPD7SS5J8q9JntDO3we4bmC99e28e0lycpI1SdaMjY0NuVxJmr+G1vWU5EvAnpMsekO73d2AQ4EnAOcmeTiQSdavyd6/qlYCKwFGR0cnXUeStPWGFhRVdeTmliV5MfCpqirgG0nuAhbRtCD2HVh1KXDDsGqUJE2tr66nzwBPAUjyKOABwI+B1cAJSR6YZD9gf+AbPdUoSaK/Gxd9BPhIku8BdwAnta2Ly5OcC1wBjAMv9YgnSepXL0FRVXcAz93MsrcBb9u+FUmSNsczsyVJnQwKSVKn+R0UGzbA4YfDrrs2jxs29F2RJM048zsoVqyAiy+GjRubxxUr+q5Ikmac+R0Ua9fC+HjzfHy8mZYk3cP8Dorly2GkPfBrZKSZliTdw/wOilWr4LDDYMGC5nHVqr4rkqQZp68T7maGJUvgoov6rkKSZrT53aKQJE3JoJAkdTIoJEmdDApJUieDQpLUyaCQJHVKcxuI2S3JGPCjrXiLRTQ3TpoNrHU4rHV4ZlO9863Wh1XV4qlWmhNBsbWSrKmq0b7rmA5rHQ5rHZ7ZVK+1Ts6uJ0lSJ4NCktTJoGis7LuA+8Bah8Nah2c21Wutk3CMQpLUyRaFJKmTQSFJ6jSvgyLJ0Ul+kGRdklP7rmdzkuyb5CtJrkxyeZJX9l3TVJLskOTbSf6571qmkmRhkvOSfL/9GR/Wd02bk+Qv2r+B7yX5RJId+65pQpKPJLk5yfcG5u2e5IIkV7WPu/VZ46DN1Pvu9u/gsiSfTrKwzxonTFbrwLLXJqkki4a1/XkbFEl2AP4WeBpwAPDsJAf0W9VmjQOvqarHAocCL53BtU54JXBl30VM0/uAf6mqxwC/ywytO8k+wCuA0ao6ENgBOKHfqu7hY8DRm8w7FbiwqvYHLmynZ4qPce96LwAOrKqDgB8Cp23vojbjY9y7VpLsCxwFXDvMjc/boAAOAdZV1dVVdQdwDnB8zzVNqqpurKpvtc830nyR7dNvVZuXZCnwdOD0vmuZSpJdgcOBMwCq6o6quqXfqjqNAA9KMgLsBNzQcz2/VVUXAT/dZPbxwJnt8zOBZ27XojpMVm9VfbGqxtvJrwNLt3thk9jMzxbgvcDrgKEelTSfg2If4LqB6fXM4C/fCUmWAY8HLum3kk5/Q/PHe1ffhUzDw4Ex4KNtV9npSXbuu6jJVNX1wF/T7D3eCNxaVV/st6opLamqG6HZ4QH26Lme++IFwOf7LmJzkhwHXF9V3xn2tuZzUGSSeTP6WOEkuwCfBF5VVbf1Xc9kkhwL3FxVl/ZdyzSNAAcDH6yqxwO/YGZ1j/xW279/PLAfsDewc5Ln9lvV3JTkDTRdvmf3XctkkuwEvAF44/bY3nwOivXAvgPTS5lBzfhNJbk/TUicXVWf6rueDk8CjkvynzTdeU9J8vF+S+q0HlhfVRMttPNogmMmOhK4pqrGquo3wKeAJ/Zc01Q2JNkLoH28ued6ppTkJOBY4Dk1c080ewTNDsN32v9rS4FvJdlzGBubz0HxTWD/JPsleQDNoODqnmuaVJLQ9KFfWVXv6bueLlV1WlUtraplND/TL1fVjN3rraqbgOuSPLqddQRwRY8ldbkWODTJTu3fxBHM0IH3AauBk9rnJwHn91jLlJIcDZwCHFdVv+y7ns2pqu9W1R5Vtaz9v7YeOLj9e97m5m1QtANWLwO+QPOf7dyqurzfqjbrScDzaPbO17b/jum7qDnk5cDZSS4DlgN/1XM9k2pbPecB3wK+S/P/d8ZcciLJJ4CLgUcnWZ/khcA7gKOSXEVzdM47+qxx0Gbq/QCwALig/X/2970W2dpMrdtv+zO3ZSVJmgnmbYtCkjQ9BoUkqZNBIUnqZFBIkjoZFJKkTgaFJKmTQSFJ6mRQSEOQ5AntPQ12TLJzew+JA/uuS9oSnnAnDUmSvwR2BB5Ecz2pt/dckrRFDAppSNpriH0T+DXwxKq6s+eSpC1i15M0PLsDu9BcO2jG3LJUuq9sUUhDkmQ1zaXW9wP2qqqX9VyStEVG+i5AmouSnAiMV9U/tvdn/1qSp1TVl/uuTbqvbFFIkjo5RiFJ6mRQSJI6GRSSpE4GhSSpk0EhSepkUEiSOhkUkqRO/x84Bz2OkEDgTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "means = [[-57.51107027  32.98489643]]\n",
      "covariances = [[ 90.24987882 429.45764867]]\n",
      "weights =  [[0.13317238 0.86682762]]\n"
     ]
    }
   ],
   "source": [
    "for i in range(0,len(labels)):\n",
    "    if labels[i] == 0:\n",
    "        plt.scatter(i, data.take(i), s=15, c='red')\n",
    "    elif labels[i] == 1:\n",
    "        plt.scatter(i, data.take(i), s=15, c='blue')\n",
    "plt.title('Gaussian Mixture Model')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()\n",
    "print(\"means =\", gmmModel.means_.reshape(1, -1))\n",
    "print(\"covariances =\", gmmModel.covariances_.reshape(1, -1))\n",
    "print(\"weights = \", gmmModel.weights_.reshape(1, -1))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
